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Abstract 

Motivated by recent observations in neuronal systems we investigate all-to-all networks of non- 
identical oscillators with adaptive coupling. The adaptation models spike-timing-dependent plas- 
ticity in which the sum of the weights of all incoming links is conserved. We find multiple phase- 
locked states that fall into two classes: near-synchronized states and splay states. Among the near- 
synchronized states are states that oscillate with a frequency that depends only very weakly on the 
coupling strength and is essentially given by the frequency of one of the oscillators, which is, however, 
neither the fastest nor the slowest oscillator. In sufficiently large networks the adaptive coupling is 
found to develop effective network topologies dominated by one or two loops. This results in a multi- 
tude of stable splay states, which differ in their firing sequences. With increasing coupling strength 
their frequency increases linearly and the oscillators become less synchronized. The essential fea- 
tures of the two classes of states are captured analytically in perturbation analyses of the extended 
Kuramoto model used in the simulations. 



I. INTRODUCTION 



Understanding the collective dynamics of coupled oscillators is an important issue in non- 
linear dynamics. In particular the coherence and synchronization of oscillators is relevant in 
many areas of science and technology. Well-studied physical examples are arrays of Joseph- 
son junctions (e.g. 113 2(1 ) and lasers (e.g. Hi), where synchronization is often desired since it 



can enhance the output power of devices. The understanding of synchronization of oscillators 



has also informed the development of control for groups of self-propelled agents 12911 . Clas- 



sical biological examples for oscillator arrays are networks of neurons. The coherence and 
synchronization of neural spiking within such ensembles of neurons underlies various types 
of rhythmic activity, which have been associated with a variety of brain functions Thus, 
the communication between different brain areas can be enhanced during certain phases of 
their rhythms, which may allow to limit effective communication to areas whose rhythms are 
(transiently) coherent 11311 . Rhythms like theta- or gamma-oscillations can provide a 'clock' 
that allows information to be encoded in terms of the timing of neuronal spikes relative to the 
phase of the ensemble oscillation |2^]. For various brain functions it has been reported that 
the relevant information is carried by correlations between neuronal spiking rather than by 



their mean firing rate la |33H . In other situations it is not desired that neurons fire in near 
synchrony but rather in a specific sequence. This is, for instance, the case for networks serv- 
in g as central pattern generators that control the movement of limbs in legged locomotion 



Il6l. ll4f] and for networks controlling the production of bird songs 111 111 . 

A unified description of the dynamics of coupled oscillators is possible if the coupling is 
sufficiently weak. The interaction between oscillators affects then predominantly their phase 



and the system can be described as a network of phase oscillators IllOl. 



tion is determined by the phase resetting curve |l 



0, 



Their interac- 



1611 . which results from the impact of the 



synaptic coupling on the dynamics of the individual oscillators. In the limit of weak coupling 
the interaction simplifies significantly and depends only on the difference between the oscil- 



3111 . in which 



lator phases. A minimal model of this type is the classic Kuramoto model H20L 
the interaction is taken to be purely sinusoidal. It has provided an excellent framework for 
understanding the onset of synchronization in globally coupled networks of oscillators with 
different natural frequencies. 

In particular in biological systems the properties of the interacting elements themselves 
as well as their interactions need not be constant in time; often they evolve on slower time 
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scales in response to the dynamics of the system. In networks of neurons synaptic plasticity, 
i.e. the modification of their coupling strengths, is a widely observed mechanism that endows 
the system with the ability to learn, to memorize, and to adapt to variable environments. 
In one well-established type of synaptic plasticity the modification of the coupling strength 
depends on the timing of the pre-synaptic input and the post-synaptic activity. Typically, the 
coupling strength is potentiated if the pre-synaptic neuron provides synaptic input before the 
post-synaptic neuron spikes, while in the converse case the synaptic strength is depressed 
For neural oscillators this tends to enhance the impact of faster oscillators on the slower ones 
and weaken the converse influence. The effect of such a spike-timing dependent plasticity 
(STDP) on the synchronization of (neural) oscillators has been studied by a number of authors 
employing extensions of the Kuramoto model. It was found that the plasticity can enhance 



the synchronization of oscillators 111811 . Moreover, for coupling strengths that are sufficient to 
render all oscillators phase-locked to each other this type of plasticity was found to lead to 
only a single completely phase-locked state. Its effective network structure has no loops and 



its frequency is given by that of the fastest oscillator 112111 . 

Synaptic plasticity need not always be homosynaptic, i.e. the modification of the strength 
of a given synapse need not depend only on the activity of the neurons connected by that 
synapse. Instead, various situations have been identified in which the strength Kij of the 
synapse from neuron j onto neuron i is modified also in response to the activity of other 
neurons I / j that synapse onto neuron i (heterosynaptic plasticity). In particular, it has 
been found that the potentiation of synapse Ku can lead to the depression or depotentiation 
of synapse Jl2]. In addition, in some preparations also the converse was observed; there 
depression of synapse K a led to the potentiation of synapse 112711 . Moreover, for that 
system evidence was presented that suggested that the sum of the weights of all incoming 



synapses remained essentially constant despite the changes in the individual synapses [27]. 
Recently, similar observations were made on an anatomical level, where the combined size of 
all synapses on a dendritic segment was found to be constant, while individual synapses grew 
or shrank in response to potentiating stimuli (3D. 

Motivated by the observation of heterosynaptic plasticity that approximately preserves 



the total weight of all incoming synapses 112711 . we investigate here a minimal model of neural 



oscillators with spike-timing dependent plasticity that conserves the total incoming weight. 
Heterosynaptic plasticity introduces competition between the synapses and the weight con- 
servation implies that even the fastest oscillator, which ends up without any inputs in the case 



3 



of the usual STDP rule, receives inputs and the network of effectively coupled oscillators de- 
velops loops. We find that this leads to qualitative changes in the dynamics of the network. In 
particular, we find not only a single state in which all oscillators are phase-locked but a host of 
such states. They fall into two classes: near-synchronous states and splay states. Depending 
on the shape of the plasticity function we find a number of different near-synchronous states, 
characterized by different dependencies of the frequency on the overall coupling strength. 
Among them are states whose frequency is essentially given by that of one of the oscillators 



in the network. In contrast to the case of purely homosynaptic plasticity 12111 this is, however, 
not the fastest oscillator but an intermediate one. In the phase-locked splay states the phases 
are distributed quite homogeneously over the whole range [0, 2ir]. While typically the order 
parameter that characterizes the synchronization of the oscillators increases with coupling 
strength, in these splay states it decreases and the oscillators become less synchronized with 
increasing coupling strength. In a neural context splay states correspond to states of the 
network in which the neurons fire in sequence spread over the whole period of the network 
oscillation. We find that the firing sequence of the splay states depends sensitively on the 
initial conditions, leading to a large number of stably coexisting splay states differing in their 
firing sequence. 



The splay states exhibit parallels to the states with sequential firing obtained in llllll . 
There it was found that networks of excitable neurons with a related type of heterosynaptic 
plasticity can produce firing sequences that match important aspects of the neural activity 
observed during the production of bird songs. 

The paper is organized as follows. In Section 2 we discuss the oscillator model and its 
connection to general oscillators and we introduce a plasticity rule that reflects spike-timing- 
dependent plasticity and conservation of incoming weights. In Section 3 we consider networks 
with few oscillators and complement the numerical simulations with a perturbation analysis 
that reveals the origin of the transitions between different phase-locked regimes. In Section 
4 we investigate larger networks. They allow a multitude of different phase-locked states, 
including many stably coexisting splay states. We capture the characteristics of the simplest 
splay states with another analytical perturbation calculation. Conclusions are presented in 
Section 5. 
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II. THE MODEL 



We consider a network of N oscillators with plastic interactions in which the sum of all 
incoming weights is conserved. For the form of the interaction we assume that for sufficiently 
small frequency differences pairs of oscillators phase-lock close to synchrony. For weak cou- 
pling such a network can be described in terms of the phases 0, of the oscillators, 

1 - 

9i = Ui ~ n ^ KijHij ^ ~ ^ ' i = l,---,N, (1) 
with the 27r-periodic interaction function H(9{ — 9j) depending only on the phase differences 



lid 



1611 . In the following we assume u)i < ojj for i < j. While we allow the oscillators to have 
different natural frequencies we assume for simplicity that they are identical in all of their 
other properties. In particular, we assume that they all have the same interaction function, 
Hij(A9) = H(A9) and the same value of the sum of all incoming weights, 

N 

K= £ Ktj. (2) 

The focus of this paper are solutions in which the oscillators are phase-locked to each other 
with small phase differences. The existence and linear stability of those states is affected 
only by the leading-order expansion of H{9) around 9 = 0, H(A9) = h^ ' + h^A9 + h.o.t. Since 
the sum of all incoming weights is conserved, the contribution can be absorbed in the 
frequency of each oscillator, — > cj{ — h^°'K. Since pairs of oscillators are assumed to phase- 
lock close to synchrony for small frequency differences the coefficient has to be positive 
and can be absorbed into Ka. As a minimal model for the phase evolution we therefore use 
the classic Kuramoto-model [2( 



3 111 , which has the same leading-order behavior in 



3> 



9i = oJi-^ J2 K i:j sin (9i - 9 j), i = l,...,N. (3) 

i+j=\ 

Even for systems with general interaction functions H(A9) the Kuramoto model will capture 
the existence and linear stability of solutions in which all phase differences are small. Their 
basins of attraction will not be properly represented, however, nor will be solutions that are 
characterized by 0(1) phase differences. 

For the modifiable interactions we consider coupling strengths that evolve depending 
on the phases of the interacting oscillators, 

rk, = f (K ljA , 0,) - Ki P^ {K l 9lA) . (4) 



The weight evolution of a single synaptic connection would be given by f (Kij,6i,6j). The 
second term in © expresses the conservation of the total weight of all incoming connections 
of an oscillator. Instead of this instantaneous conservation one could also consider achieving 



homeostasis of the total weight on a longer time scale [35]. The existence of the phase-locked 
states that we are interested in here would not be affected by such a slower evolution, since 
they correspond to fixed-points of ©. At most, such a delayed homeostasis could influence 
their stability 

We assume that the weights evolve on a slow time scale, r » 1, and change only little 
during one period of oscillation of the interacting oscillators. Due to averaging the weight 



changes depend then to leading order only on the phase difference Oi — Oj 11511 . For the plastic- 
ity function f{Kij,9i — Oj) we use a functional form that is motivated by the widely observed 
spike-timing dependent plasticity (STDP) of neuronal oscillators (2D. There a synaptic con- 
nection is potentiated when the pre-synaptic neuron spikes before the post-synaptic one and 
depressed otherwise. Within the present phase framework this corresponds to a potentiation 
when the phase of the pre-synaptic oscillator is larger than that of the post-synaptic oscillator. 
Specifically we use 



f (Kij,9i,6j) 



(a-Kij)e t p 
/3 + p 1 fa - 9j 



for Oi - Oj £ (-7T, -tp) 
for Oi - Oj G [-tp, tp] , 
for Oi - Oj £ (tp, 7r] 



(5) 



where Oi — Oj is taken modulo 2tt within the range (— ix, it]. We include a central phase window 
[—ip, tp] within which potentiation and depression are continuously interpolated OJ]. Typically, 
we will assume this window to be narrow, tp <C r dp or even tp = 0. The coefficients /?o,i are 
given by 

1 



A) 



l _± 

-e (a - Kij) 



-Kije T d, 



1 

2i) 



i.i 



a) e t p 



K^e 
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Figure 1: Plasticity function f(Kij,6i - 9j) for a = 10, ip = 0.02 , r d . p = 0.1, = 5.89. The coupling 
strength is increased if the phase 0j of the presynaptic oscillator is larger than the phase 6* ; of the 
postsynaptic oscillator, A9 = Q l - 6j < 0. 

Our main control parameter is the sum K of all incoming weights. The parameter a sets 
the maximal strength of an individual synapse in the absence of the homeostatic term in ([5]). 
We focus here on phenomena that are dominated by the limitation of the overall coupling K 
and choose a well above K. Note, however, that due to the homeostatic, second term in (|4]> 
the coupling strengths Ky are not strictly limited to Kij < a. Correspondingly, we did not 
find qualitatively different behavior when a was chosen somewhat below K. 



III. FEW OSCILLATORS 



For oscillator networks in which the plastic coupling is not conserved it was found that for 
arbitrary network sizes there is only a single phase-locked state and the transition from the 
incoherent states to that phase-locked state i s hy steretic only if the plasticity windows t p ^ 
for potentiation and depression are not equal [21]. We find that with the conservation of the 
overall input strength hysteresis arises even with equal plasticity windows. Moreover, the 
transition scenario and the extent of hysteresis depends strongly on the natural frequencies 
of the oscillators. The case of three oscillators is illustrated in FigHJ Depending on ui2 with 
wi 5 3 fixed, all three oscillators can phase lock in what seems a single hysteretic transition 
(1.2 < oj 2 < 1.5) or in two subsequent transitions with an intermediate, partially phase-locked 
state (1.6 < u 2 < 1.9). 
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Figure 2: Transition sequences to the phase-locked states for N = 3 oscillators for different values of 
the intermediate frequency uj 2 - Parameters: r = 20, r d = t p = 0.3, wi = 1, L03 = 2, a = 100, iJj = 0. 
For uj 2 > 1.5 the frequency of the phase-locked state is very close to w 2 (cf. FigJS). Red symbols denote 
increasing K, black symbols decreasing K. 

A particularly striking aspect of the simulations shown in Figj2] is that the frequency uj of 
the completely phase-locked state exhibits two different regimes: for uj2 closer to the lower fre- 
quency u>\ the frequency u depends only little on u) 2 , while for oj 2 closer to the larger frequency 
W3 the three oscillators oscillate at a frequency that is very close to the natural frequency uj 2 
of the second fastest oscillator. This is shown more explicitly in FigJH We find this surpris- 
ing selection of the frequency of the second-fastest oscillator also in simulations with more 
oscillators (see below). 
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Figure 3: Dependence of the frequency of the phase-locked state on the frequency of the second oscil- 
lator, cj 2 , for K = 3, r = 20, r d = t p = 0.3, w a = 1, oj 3 = 2, a = 100, 7/) = 0.005 (cf. Fig©. The dashed and 
the dotted line give the analytical results dl4H and dl9l > for the frequency of PL0 and PL2, respectively 

To get an analytic understanding of the regimes found in FigJUand in particular to identify 
the origin for the phase-locking at a frequency close to that of the second-fastest oscillator we 
consider the situation in which the phase differences AOij are sufficiently small to allow a 
linearization of the nonlinearities in ([T]> and ©. We therefore assume that the differences 
between the three frequencies are small, 

uj 2 = oj\ + ef2 2 , ^3 = ^1 + 6^3, £«1, (6) 

so that a coupling of 0(1) can lock the phase differences at small values. The phase differences 
can therefore be expanded as 

A9 12 = e 50$ + e 2 60$ + 0(e 3 ), A8 23 = e 59$ + e 2 59$ + 0(e 3 ). (7) 

In addition, to avoid that all phase differences fall in the central range [—ip, ip] of the plasticity 
function we assume that range to be narrow, 

V> = e*. (8) 

We also expand the coupling coefficients, 

Ka = Kf + eK$ + 0(e 2 ) (9) 

The piecewise definition of the plasticity function f (Kij,9i,0j) in © requires that one 
distinguishes different cases depending on A9ij = 9i — 9j. For S7 3 — fl 2 = (003 - o; 2 )/e and 



^2 — ^i = (^2 — not too small both phase differences A#i2 and A6>23 fall outside the inner 
range [—tjj, ip] of the plasticity function. Inserting the expansions ( I7I9D into d3|4D leads then in 
a straightforward fashion to evolution equations for the leading-order contributions and 



n 2 + iwg { J$> - 2K} + '-SOS {*g> + - k} , (10) 

o 2 + lae® {k - - *g>} - iw« {irff + ^} (ii) 



and 



tK {0) 
TIX 32 



tK® = a {K-2K[f}. 



and an additional equation for one of the phases, 9\ say, which for steady states yields the 

L32 is slower than that of 



oscillation frequency oj = 0\. Note that the evolution of K$ is slower than that of KLy and 



k[^. These equations have two fixed points, which represent phase-locked solutions. Only 
one of them is linearly stable. It is given by 

PLO : K$ = K[f = \K K$ = (12) 

60$ = eus^h > % 6e W = snz+p, > ^ (13) 

with the remaining coupling coefficients determined through the conservation law, e.g. = 
K — K%i ■ This solution is only valid as long as 66$ > * and 59^ > ^ as indicated by the 
inequalities in ( fT3l ). The oscillation frequency of PLO is given by 

u = u l + e{hl- i + hl 2 ^j+0{€ 2 ). (14) 

This frequency is not necessarily close to U2\ in fact, it varies only weakly with u; 2 - 

When $7 3 — 2^2 < the phase difference A9$ of PLO falls into the central range 

[—tp, V>]. This modifies the evolution equations for (cf. ® and ( IA1IA2IA31 ) in the appendix) 
and the expansion 07191 ) yields 3 possible phase-locked solutions. They are given by 

PL1: K$ = K K[f = \K K$ = (15) 

xaW — 3 03—02 iT/ AflW — 3 \ .t, 

°°32 — 2 k — 21 — 4 k - 
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with frequency 



(16) 



and 



PL2a: 



< K$ = k 




(17) 



(18) 



with frequency 




(19) 



and a solution PL2b that is obtained from PL2a by interchanging oscillators 2 and 3, keeping 
in mind that interchanging O2 ^3 implies ^ = K — . Again, the range of 
validity of each solution is indicated by the various inequalities. 

The ranges of validity of the solutions PLO, PL1, and PL2a,b are mutually exclusive. In 
particular, depending on the sign of O3 — at most one of PL2a and PL2b is valid. Moreover, 
at the validity limit of the solution PLO, SI 3 — 2O2 = 5K^/6, it becomes equal to PL2a with 
= 0, which at the same time represents one limit of validity of PL2a. At the other limit of 
validity of PL2a one has = K. There it coincides with PL1 at one of its limits of validity. 
Finally, PL1 reaches its other limit of validity when 59 2 \ = To continue the solutions into 
this regime a further expansion would be necessary, in which also 59$ is assumed to be in 
[— Thus, we find a single branch of near-synchronous phase-locked solutions, which 
exhibit, however, quite different behaviors in the different regimes. 

Note, that in none of the regimes the oscillators are truly synchronous, i.e. their phase 
differences do not vanish. This is to be contrasted with previous results on oscillator network 
models with homo-synaptic plasticity where it had been found that the plasticity can lead 
to perfect synchronization of the oscillators, although the oscillators have different natural 
frequencies I118ll . In that model the plasticity can effectively induce different values of H(Q) 
for the different oscillators, which can compensate for the differences in natural frequencies 
even for A9ij = 0. With conserved total incoming weights, however, the plastic modification 
of H(0) is the same for all oscillators and perfect synchrony cannot be achieved. 

The quantitative comparison of the perturbation analysis and the numerical simulations 
presented in Fig|3] shows that PLO and PL2a capture the phase-locked states obtained in 
FigHJ Thus, the perturbation analysis reveals that phase-locking of the three oscillators at a 
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frequency very close to that of the second fastest oscillator is obtained if the transition region 
between potentiation and depression is narrow, ip -C 1. 

To investigate the additional transition from PL2a,b to PL1 that is predicted by the per- 
turbation calculation we perform numerical simulations that include a significant central 
range of the plasticity function, i.e. V 0. To allow a quantitative comparison with the per- 
turbation expansion we use small frequency differences. The resulting coupling coefficients, 
phase differences, and frequencies, are shown in FigJH as a function of K. The solutions are 
most easily identified by their coupling coefficients and phase differences. For small K one 
finds PLO, which is characterized by -K32 = K21 = and A#2i, A032 > V- As K is increased 
the phase difference A#32 decreases and eventually falls into the range [— ip,ip], marked by 
a dotted line in FigHJj. At this point the solution PLO transforms into PL2a and ^32 starts 
to deviate from 0. Since t\> 7^ the frequency of PL2a is not independent of K in contrast to 
what was found in the simulations shown in FigHk. In fact, relative to the small frequency 
differences used here the independence of the frequency is quite pronounced. 

As K is increased further K32 reaches the value K. There PL2a transforms into PL1. For 
yet larger values of K Figj4] reveals an additional continuous transition to a state PL3 in 
which also A#2i enters the region [—ip, ip] and K21 becomes non-zero. We have not performed 
the additional modification of the expansion to capture this state analytically. 
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Figure 4: Continuous transitions of the phase-locked state for N = 3 oscillators. Analytical results 
are denoted by dashed lines, a) Coupling coefficients K l3 as a function of K. Dotted vertical lines 
indicate the transitions between different regimes, b) Phase differences. The border of the central 
region [— ip, tp] of the plasticity function is indicated by a dotted line, c) Frequency of the phase-locked 
state as a function of the overall coupling K. Parameters: a = 10 (thick lines) and a = 20 (thin solid 
lines), r = 100, wi = 1, ui 2 = 1.03, u> 3 = 1.1, r d , p = 0.1, ip = 0.02. 



The frequency and coupling coefficients obtained from the perturbation calculation 
(dashed lines in Fig© agree quite well with the numerical simulations (thick solid lines). 
Nevertheless, for PL2a the differences are quite noticeable. While the upper limit a of the 
individual synaptic strengths does not appear in the leading-order results ( I17I19D of the per- 
turbation calculation, it turns out that contributions proportional to a" 1 arise at next order, 
which become large for small a (cf. Appendix |A|>. Thus, increasing a from a = 10 to a = 20 
further improves the agreement (thin solid lines in Figj4]). 

Thus, even in this system comprised of only 3 oscillators the combination of the central 
range [— ip, ip] of the plasticity function with the conservation of incoming coupling strengths 
leads to transitions between at least four regimes in which the phase-locked solution exhibits 
quite different behavior. As noted before, the transitions between these regimes do not repre- 
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sent bifurcations associated with instabilities but points at which the plasticity function © 
in the underlying differential equations is not differentiable or a coupling strength reaches 
the maximal value imposed by the weight conservation. 

IV. MANY OSCILLATORS 

For larger networks of oscillators an additional, qualitatively different class of stable 
phase-locked states arises. Sample transition sequences for increasing and for decreasing K 
are shown in Figj5]for N = 20 oscillators with frequencies equally spaced in the interval [1,2]. 
In both cases we start with homogeneous coupling, Kij = j^zfK, but random phases. For in- 
creasing K the initial spread in the frequency of the unsynchronized oscillators decreases and 
step by step the six fastest oscillators merge into a cluster oscillating with a single frequency. 
At K = 29 all oscillators phase-lock and form a new state, the frequency of which is higher 
than the natural frequency of the fastest oscillator and increases further with increasing cou- 
pling. This new state persists to the largest values of K investigated. Decreasing K from 
large values - again starting with homogeneous coupling - a different globally phase-locked 
state is reached. Its frequency is very close to that of the third fastest oscillator. Near K = 50 
it crosses over to a phase-locked state with a frequency very close to that of the second fastest 
oscillator. At K = 40 that state undergoes a jump transition to the phase-locked state found 
when increasing K from small values. 



14 




2 3 
Phases 



Figure 5: Transitions between three different, globally phase-locked states for N = 20 oscillators, a) 
Frequencies as a function of K. There are two PL2-like states with frequencies close to those of the 
second- and third-fastest oscillator, respectively, marked by dashed-dotted lines. The analytical result 
d29l > for the splay state is marked by a dashed line, b) Matrices of the coupling coefficients for the two 
PL2-like states (bottom panels, K = 59.4 and K = 45.5) and the splay state (top panel, K = 36.4). c) 
Order parameter rasa function of K (cf. H20D). It increases with increasing K for PL2-like states, but 
decreases for the splay state, d) Phases of the two near-synchronous PL2-like states and of the splay 
state. Parameters r = 20, t p = 0.15, Td = 0.3, a = 100, ip = 0. 



To understand the main aspects of the phase-locked states consider the coupling coeffi- 
cients Kij established by the plasticity. With only homosynaptic plasticity each oscillator 
would be coupled with the maximal strength a to all faster oscillators and would receive 
no input from any of the slower oscillators. The magnitude of the phase difference between 
the oscillators would affect only how fast these final values of the coupling coefficients are 
reached. The heterosynaptic plasticity employed here introduces competition between the in- 
coming couplings and the resulting steady-state values are distributed over the whole range 
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[0,-K]. If t p is small the input to each oscillator is predominantly coming from a single other 
oscillator (Figj5]b). This allows to define chains of dominant coupling. Since the conserva- 
tion of the overall incoming coupling enforces that each oscillator receives input, these chains 
must contain loops. 

The coexisting phase-locked states differ in characteristic ways in their chains of dominant 
coupling. Similar to the state PL2 given by 017I19D the states oscillating with frequencies close 
to those of oscillators 2 and 3, respectively, have strong input from oscillator 2 to oscillator 1 
(Figj5p bottom panels). They differ in additional input from oscillator 3 into oscillators 1 and 
2. Although these two coupling coefficients are relatively small, they are sufficient to pull the 
frequency down to that of oscillator 3. They go to in the cross-over near K = 50. In both 
states all remaining oscillators receive their dominant input from a single other oscillator, 
which in most cases is the oscillator with the next higher frequency. Thus, in both states the 
chain of dominant coupling contains only a small loop involving the fastest oscillators 19 and 
20 or 18, 19, and 20, respectively. In the jump transition at K = 40 the input from oscillator 
2 into oscillator 1 disappears and instead oscillator 1 receives strong input from the second 
slowest oscillator (Figj5]b top panel). Consequently, in this state the chain of coupling consists 
of a single large loop involving essentially all oscillators. 

The qualitative difference between the different types of phase-locked states manifests 
itself also in their phase differences. While in the PL2-like states the phases are closely 



clustered, the phases of the other state are distributed 



32, 



quite 



34, 



lomogeneously over the interval 



[0, 2tt] (FigUJi) identifying it as a splay state [30|, 
differ therefore significantly in terms of the order parameter 

N 



36]. The states in the two classes 



1 

N 



(20) 



which characterizes the degree of synchronization of the state. Typically one would expect 
that the synchronization becomes stronger as the coupling between the oscillators is in- 
creased. This is indeed the case for the PL2-type states (FigHh). However, in the splay state 
the order parameter decreases with increasing coupling, indicating that the coupling tends 
to spread out the phases more uniformly. 



16 



A. Perturbation Analysis 

To understand the origin of the splay state we again employ a perturbation analysis. It 
is guided by the observations shown in Figj5l The characteristic features of the splay state 
can be captured by considering a regime in which each oscillator interacts only with one other 
oscillator. This is the case if the window for potentiation t p is sufficiently small and the phases 
are distributed sufficiently homogeneously. Thus, we assume 

t p < min (Afli+1,0 (21) 

' l<i<N ' 

and 



maxA0 i+ i j < minA6>j +2 ,i. (22) 

i i 

To allow the linearization of the equation of motion for the phases we assume in addition that 
the number of oscillators is large, 

N > 1, (23) 

so that A0 i+M = 0{N- 1 ). 

This perturbation analysis will be strictly valid in the limit N — > oo, which implies that 
all phase differences lie in the central region [—ip, ip] of the plasticity function. We expect, 
however, that this approach will also give good results for intermediate values of N for which 
minj A^j+i^ > ip. For simplicity we therefore take in this analysis ip = with the expectation 
that the results will also apply to systems with tp > as long as N is not too large and 
therefore miiij A0j + i j > ij). 

Here and in the following the phase indices are considered modulo N. Thus, in particular, 
A(9jv,./v+i = A0jvi- For the splay state of interest we assume that A6m - 2vr = On — #i - 2-7T = 
0(l/N). 

Independent of the assumption d2H ). \A6ij\ > ip implies that for j < i eq.® always has a 
solution Kij = 0. For j > i the equations © for Kij are simplified by the assumptions d21|22D . 
which imply 

_ 1 A A _lAA 

e l+m ' l < e t p f orm >2. (24) 
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Thus, 



A 



1.1+ in 



(a - e-p + hr— < - 2^ K M e d +Z^( a - K W e p X25) 



A' 

(a - A M+m ) e i A ^+™ - |( a _ ei A ^ +1 + fc.o.i.} . (26) 



For m = 1 the two terms are of the same order and with a > A one obtains for the fixed point 
Ki i+ i = A. For m > 2 the first term can be neglected relative to the second one due to ( 1241 ) 
and one has Aj j +m = 0. In summary, to leading order the coupling coefficients for the steady 
state are given by 

A M+ i = A i = l,...,N-l, 
K N1 = A, (27) 
tfij = j^i + l,i^N. 

For the phase differences one obtains from dU) for the phase-locked state oscillating with 
frequency u 

A0 M+1 = 2 (Wi - w) (28) 
A 

to leading order in A^+i. This direct connection between the phases and the frequencies 
shows that condition d22l amounts to the assumption that the natural frequencies are not 
distributed too heterogeneously. 

The common frequency u of the oscillators is obtained by expressing A0/vi in two ways. On 
the one hand one has 

N-l N-l 

A6 N1 = 0n-O 1 = -V] A0 M+i = —^Y\ (wj - w) . 

i=i K i=i 

On the other hand, using A8 m - 2vr = 0(1/N) in © with i = N yields 

iV 

A6*tvi = 2ir + — (wat — w) . 
A 

Combining the two expressions for A9ni results in 

K . , I 

uj = oj + 2tt— tt with cj = — > cji. (29) 
iv^ N z — ' 

i=l 



Replacing w in d28l ) the phase differences are given by 



A0 M+1 = -^ + ^ (30) 
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Our analysis assumed A^iat > 0. With 028I29D this implies that the splay state exists only 
above a minimal coupling strength K c , 

N 2 

K>K c = — (oj n -uj), (31) 

and its frequency is above that of the fastest oscillator, uj > uj^. 

Within the framework of ( [24] ) small perturbations to the coupling coefficients Kij decouple 
from the perturbations of the phases 9i and it is easy to show that the splay state is linearly 
stable as long as a > K. 

The analytical result < f30b shows that with increasing K the phases become more evenly 
distributed, independent of the natural frequencies of the oscillators. This results in a de- 
crease of the order parameter r with increasing K in agreement with the numerical simula- 
tions (Figj5b)- Eq.d29l captures the linear growth of the oscillation frequency with K (dashed 
lines in Figj5k). For the parameters of FigHk the agreement is, however, not quantitative. 
In the analytical calculation we considered the limit of small t p , which allows to assume that 
each oscillator receives inputs only from a single other oscillator. In Figj5]3 this is not quite 
the case. Reducing the plasticity window to t p = 0.05 with = 0.1 yields, however, very 
good quantitative agreement (Fig©. Again we find extensive bistability between the splay 
state and a PL2-like state oscillating with the frequency of the second-fastest oscillator (dash- 
dotted line). For these parameters we found no transition from the PL2-like state to the splay 
state when decreasing K. 




A 

Overall Coupling K 



Figure 6: Quantitative agreement of analytical and numerical results for the splay state for N = 25 
oscillators with t p = 0.05 and r d = 0.1. Analytical result ([29} denoted by dashed line. Frequency of 
second-fastest oscillator denoted by dashed-dotted line. Other parameters: r = 20, a = 500, ijj = 0. 
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The fact that the oscillation frequency of the splay state is larger than the natural fre- 
quency of the fastest oscillator can be seen to be a direct consequence of the conservation of 
total incoming weights. It induces an input from the slowest to the fastest oscillator. Since 
for sufficiently large N the slowest oscillator lags the fastest one by almost 2ir, the slowest 
oscillator is effectively pulling the fastest oscillator ahead. 

B. Multiplicity of Attractors 

Figs J5l6l show extensive bistability between splay states and PL2-like states. Moreover, 
Figj5p,d shows that this splay state does not exactly correspond to the analytically obtained 
solution since the coupling sequence of oscillators 1 to 6 and with it their firing sequence does 
not strictly follow their natural frequencies. This suggests that splay states with other firing 
sequences may exist stably as well. 

To investigate the multiplicity of attractors for these larger oscillator networks we have 
performed simulations with 500 different initial conditions for the phases of the oscillators, 
keeping the initial coupling coefficients homogeneous, and with different ramping rates for 
the overall coupling K. The latter is motivated by the observation that in Figs J5l6l the splay 
states were obtained by ramping K up from small values, while the PL2-like states arose 
when K was set instantly to a value in the phase-locked regime. As expected, the fraction of 
initial conditions that lead to splay states rather than PL2-like states increases with decreas- 
ing ramping rate for K (Figl7]>. 




0.01 0.1 

A 

Ramping Rate dK/dt 

Figure 7: The fraction of initial conditions leading to splay states rather than PL2-like states decreases 
with the ramping rate dK/dt. Parameters: N = 20, K initial = 30, Kf ina i = 60, t p = 0.05, r d = 0.1, 
t = 20, gl = 500, V = 0. 
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The splay states reached from the different initial conditions are not all the same. In fact, 
none of the 267 splay states obtained for dK/dt = 0.006 had the same firing sequence. This is 
apparent in the firing matrix F shown in FigJSp. where the color of the element Fy indicates 
for run i the oscillator that fired at the j th -position in the firing sequence. The rows are 
ordered by the number of the oscillator that fires first, second, third, etc. Analogously, none of 
the firing sequences of the 233 initial conditions (out of 500) that led to PL2-states appeared 
twice (FigHb). 




Figure 8: Firing sequences of the states reached using 500 random initial conditions for the phases 0, 
but homogeneous values . Each row gives the firing sequence for one initial condition, a) 267 splay 
states, b) 233 PL2-like states. Parameters as in FigPTl ramping rate dK/dt = 0.006. 

Despite their different firing sequences, the various PL2-like states oscillate with a fre- 
quency that is extremely close to that of the second-fastest oscillator. This is not the case for 
the splay states: their frequencies are quite broadly distributed (Fig®. The difference be- 
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tween the two types of state can be understood intuitively. The PL2 states are dominated by 
the fastest 2 or 3 oscillators. Different firing sequences of the slower oscillators have therefore 
little impact on the overall state. In the splay states, however, the fastest oscillator is pulled 
ahead by one of the slow oscillators, which in turn is pulled ahead by another oscillator and 
so on until the circle closes with the fastest oscillator pulling a slower one. Since (almost) 
all oscillators are part of this chain of dominant coupling, the overall state and its frequency 
depend significantly on the firing order of the slower oscillators. 




Frequency co 



Figure 9: Frequency distribution of splay and PL2-states. The splay states have a broad frequency dis- 
tribution while the frequencies of the PL2-states are indistinguishable from that of the second fastest 
oscillator (marked by a triangle). The corresponding firing sequences are shown in FigJHJ Parameters 
as in Figjjl ramping rate dK/dt = 0.006. 

Can the chain of dominant coupling contain more than a single loop? FigflO] shows that 
this is indeed the case. FigflQa depicts the coupling matrix Kij obtained for one set of initial 
conditions of the phases with ramping rate dK/dt = 0.006 (cf. Figs J8l9D in which the fastest 
oscillator O 2 o (marked by a hashed circle in FigflOb) gets significant input from two slow 
oscillators, Oq and Og (marked by solid circles in FigflOb). While oscillator Og drives only O20, 
oscillator Oq drives in addition also oscillator 0\, dividing the chain of dominant coupling 
and generating two loops. The large loop involves all oscillators and reaches eventually 0$, 
while the small loop involves only O20-18, O7, and Oq. Oscillator 8 lags O20 by almost 2tt 
and is therefore effectively pulling O 2 o ahead as in the splay state described above. Oscillator 
Oq, however, is only slightly behind O 2 o; it actually holds O20 back and increases the phase 
difference between oscillators 0$ and O20 leading to a tighter clustering of the phases of the 
fastest oscillators C>2o-i8- 
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2 3 4 
Phase 



Figure 10: Phase-locked splay state with two loops in the chain of dominant coupling, a) Coupling 
coefficients Kij. Oscillator O 2 o receives input from O e and 8 . b) Phases of the oscillators with the 
chain of dominant coupling marked by dotted lines. 6 (solid circle) couples to 0\ and O 2 o, holding O 2 o 
back. Parameters as in Fig|j] ramping rate dK/dt = 0.006. 



FigJTTl gives an overview of the splay states shown in Figj8]in terms of the strengths K 2 oj 
of the incoming links of the fastest oscillator O20 (FigJTTk) and their frequency (FigJTlb). In 
the single-loop splay states there is only a single such strong input and it has full strength K. 
In the two-loop splay states, however, at least two oscillators provide significant input to O20, 
each with smaller amplitude. With the rows in FigfTlb being sorted by increasing frequency, 
it is apparent that the bimodal structure of the frequency distribution of the splay states seen 
in Figj9] reflects the occurrence of 1-loop and 2-loop states, respectively. 
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Firing sequence position 



Frequency 



Figure 11: Phase-locked splay states with one and two loops in the chain of dominant coupling, a) 
Coupling strengths of the incoming links of oscillator O 2 o for the 267 runs resulting in splay states 
sorted by increasing frequency (cf. Fig|8}- Single-loop states have a single incoming link of strength 
~ K, while 2-loop states have multiple, weaker incoming links, b) The frequencies of these states. 
Parameters: Parameters as in FigJTl ramping rate dK/dt = 0.006. 

Thus, the synaptic competition introduced by the heterosynaptic plasticity combined with 
the weight conservation stabilizes a variety of splay states with characteristic firing se- 
quences. 

V. CONCLUSION 

We have investigated the synchronization and phase-locking of networks of weakly coupled 
oscillators whose interactions evolve in response to the dynamics of the oscillators. Specifi- 
cally, we considered coupling strengths that are modified slowly depending on the phase dif- 
ference between the oscillators involved, while keeping the total weight of the incoming con- 
nections of any given oscillator constant. This was motivated by observations in neural sys- 
tems, where spike-timing dependent plasticity is found quite commonly. Our consideration of 
heterosynaptic plasticity that conserves total incoming weight was triggered by experiments 
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in which the overall strength of all incoming synapses was found to remain approximately 
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constant while individual synapses were potentiated or depressed |3|, 

For purely homosynaptic plasticity there is only a single state in which all oscillators are 
phase-locked to each other 112 ill . In this state each oscillator is coupled equally to all faster 
oscillators and the overall frequency is that of the fastest oscillator. Including heterosynaptic 
plasticity with weight conservation, we find a host of different phase-locked states, which fall 
into two classes: near-synchronous states and splay states. 

Due to the continuous transition in the plasticity function between potentiation and de- 
pression and due to the conservation of the overall coupling the near-synchronous solutions 
exhibit quite different behaviors in different parameter regimes. This is reflected in partic- 
ular in the dependence of the oscillation frequency on the overall coupling strength. For 
the case of three coupled oscillators we identified various continuous transitions analyti- 
cally. If the transition region of the plasticity function is narrow the frequency of the near- 
synchronous solutions depends only weakly on the overall coupling strength. Interestingly, 
there are large parameter regimes in which the frequency is essentially given by that of one of 
the oscillators in the network, which is, however, neither the fastest nor the slowest oscillator. 

In the splay states the phases are distributed over the whole interval [0, 2ir]. A large num- 
ber of different stable such states are found. In simple splay states the chain of dominant 
coupling, which represents their effective network structure, forms a single loop and defines 
a firing sequence characteristic for that splay state. In addition, we also found more complex 
splay states with a 2-loop structure. A multitude of splay states with different firing se- 
quences coexist stably. Their oscillation frequencies are broadly distributed, reflecting their 
different firing sequences. Strikingly, the splay states become less synchronized when the 
coupling strength is increased. At the same time their overall oscillation frequency increases 
essentially linearly. This frequency is larger than that of the fastest oscillator: the fastest 
oscillator is pulled ahead by one of the slow oscillators. The essential aspects of the splay 
states are captured quantitatively in analytical perturbation calculations. 

The splay states are characterized by the unidirectional ring topology of their chain of 
dominant coupling. For fixed, non-plastic coupling strengths the dynamics of oscillators that 
are coupled unidirectionally in a ring has been studied in detail previously 
Such coupling leads quite naturally to oscillatory states in the form of traveling waves, which 
correspond to the splay states found here. Results for their stability have been obtained 



for the Kuramoto model and extensions thereof [26]. For unidirectionally coupled Duffing 
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oscillators their instability has been identified as an Eckhaus instability 11171 



19, 



2J. The 



effect of a delay in the interaction in such networks has also been discussed for an amplitude- 



equation model and for coupled FitzHugh-Nagumo neurons 12311 . For pulse-coupled oscillators 
the phase-locked solutions have been described using maps for the firing times 1^,0,0]. Based 
on the phase-resetting curves, these analyses showed how the stability of the phase-locked 
states can be controlled by modifying the phase-resetting curves through slight modifications 
of the neural dynamics. Thus, in networks functioning as central pattern generators the 
network dynamics can, for instance, be switched between different animal gaits by injecting a 
steady current into the neurons These results shed some light on the phase-locked splay 
states investigated here. However, while in these previous analyses the network structure 
and the coupling coefficients were kept fixed, an essential part of the dynamics discussed 
here consists of a restructuring of the chain of dominant couplings. 

We have described the network evolution in terms of the self-organization of a network of 
oscillators with different natural frequencies in the absence of any input. Once established, 
some of the splay states turn out to persist if all frequencies uji are set to the same value, 
u)i = U). Thus, eqs. dll4D can also be read as describing a network of identical neural oscillators 
that receive heterogeneous tonic input, which modifies their firing rate (natural frequency) 
and which can be used to train the network to generate different firing sequences. However, 
due to the sensitive dependence of the firing sequence on the phase distribution of the initial 
conditions it would be necessary to control also the oscillator phases during the training 
period to select specific firing sequences. 

The focus of this work was the effect of heterosynaptic plasticity on the dynamics of a net- 
work of oscillators. In particular, our model has been motivated by the experimental finding 
that in certain neurons in the amygdala heterosynaptic plasticity roughly balanced homosy- 



2711 . In other sys- 



naptic plasticity keeping the overall coupling approximately constant (3|, 
terns heterosynaptic plasticity may not conserve the overall synaptic weights. Thus, it has 
been observed that heterosynaptic plasticity can alternatively be controlled by the overall ac- 
tivity of the postsynaptic neuron, independent of its inputs 112 511 . or that it can reflect limited 
resources (proteins) of the neuron B12I1 . In particular in the latter case, it would be natural to 
model the plasticity by limiting rather than fixing the overall weight of all synapses, as has 
been done in a model for sequence generation in bird song II 1111 . 

For the description of the oscillators and their interaction we chose a phase model. This is 
adequate in the limit of weak coupling. The phase model is characterized by its interaction 
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function H(A0), which depends on details of the dynamics of the uncoupled oscillators and on 
their phase-resetting curves Thus, type-I oscillators, which arise from a saddle-node 

bifurcation on a circle, and type-II oscillators, which arise from a Hopf bifurcation, typically 
lead to different functional forms of H(A9). Similarly, type-I phase resetting curves, which do 
not change sign, and type-II phase resetting curves, which do change sign, result in different 
forms for H(A8) . All of the phase-locked states investigated here are characterized by very 
small phase differences. Therefore only the behavior of H(A9) in the immediate vicinity of 
A9 = is relevant for their existence and linear stability. Moreover, due to the conservation of 
the total incoming weights the constant contribution H(0) can be absorbed into the frequen- 
cies of the individual oscillators. The core of our results apply therefore independent of these 
different types of oscillators and phase resetting curves as long as the interaction is such that 
the oscillators phase-lock near synchrony when their frequencies are not too different, as is 
the case in the minimal form of the classic Kuramoto model. Thus, while we were mainly 
motivated by the dynamics of neural networks, addressing the issue of synchronization and 
sequential firing, we expect that our results apply to a much larger class of adaptive oscillator 
networks in the weak coupling regime. 
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Appendix A: Higher-Order Expansion for 3 Oscillators 



Here we give some more details for the expansion of the solution PL2a to order C(e 2 ), 
which reveals the dependence of that solution on a and an additional dependence on K . 

Inserting A6I7I8I9D into the 3 phase equations © results at each order in two equations for 
the phase differences S9^ 2 ' and 89^ as well as an equation for one of the individual phases, 
83 say. For phase-locked solutions the latter equation determines the overall frequency of 
oscillation via lo = #3. To wit, for PL2a one obtains at 0(e) dlOlllD and 



as well as modified equations for K 



(o) 



2^K 
1 

K 

The fixed-point solutions of these equations are given by ( fl5l ) and ( fT71 ) 



(69$ - *) , (Al) 



Tk[V = Z(k-2K®). (A3) 



To illustrate the form of the contributions at the next order we focus on PL2a. Its fixed- 
point equations read at 0(e 2 ) 
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Their solution is given by 

,(2) 
'32 

,(2) 
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This results in an overall frequency given by 

u = Ul + e (n 2 + - -e 2 — f 12fi 2 - 6O3 + 5if#) + C(e 3 ). 
V / 3 arrf V / 

Thus, the corrections at 0(e 2 ) show that the fixed-point solution is not independent of a as 
might have been assumed based on the leading-order results. Moreover, the a-dependence of 
the corrections is of 0(a~ 1 ). Thus, with increasing a the 0(e 2 )-corrections, in particular to the 
frequency, decrease, consistent with the improved agreement of the perturbation calculation 
with the numerical simulations seen in FigHl Moreover, the phase difference A# 3 2 for PL2a 
is not independent of if and does not lie exactly at the border of the central region of the 
plasticity function. 



1 if* 
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